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Abstract 

We study the quark mass dependence of the finite temperature QCD phase transition in the heavy 
quark region using an effective potential defined through the probability distribution function of the average 
plaquette. Performing a simulation of SU(3) pure gauge theory, we first confirm that the distribution 
function has two peaks indicating that the phase transition is of first order in the heavy quark limit, while 
the first order transition turns into a crossover as the quark mass decreases from infinity where the mass 
dependence of the distribution function is evaluated by the reweighting method combined with the hopping 
parameter expansion. We determine the endpoint of the first order transition region for Nf = 1, 2, 3 and 
2 + 1 cases. The quark mass dependence of the latent heat is also evaluated in the first order transition 
region. 



* We found an error in the analysis program of the effective potential developed for this paper after the publication, 
Phys. Rev. D 84, 054502 (2011). The error was in a coefficient of the term proportional to the plaquette. This 
causes slight shifts in the values of /?trans at k > and thus /3 cp . Accordingly, we have replaced FIG. [3] FIG. [9l 
TAB. HQ and Eq. |g0j from those given in the published version arXiv:1106.0974v2]. On the other hand, this error 
does not propagate to d 2 V B a / dP 2 . Therefore, the discussions and the conclusions of the paper, including the values 
of Kcp as well as other figures and tables, are not affected. 

t Present address: Physics Department, Bookhaven National Laboratory, Upton, New York 11973, USA. 
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I. INTRODUCTION 



The QCD phase transition is one of the important ingredients in understanding the evolution of 
the early universe. Monte-Carlo simulation of lattice QCD is the most powerful approach to study 
the nature of the QCD phase transition at present It is known that the order of the phase 
transition depends on the values of quark masses: The deconfinement phase transition is of first 
order when masses of all three (up, down and strange) quarks are either sufficiently large or small, 
while it turns into a crossover in the intermediate region 0-0] • Previous studies with staggered 
quarks strongly suggest that the transition becomes a crossover at physical quark masses [8Hl0l]. 
A confirmation of this result by other fermion formulations such as Wilson-type fermions or by 
other methods of analyses is, however, mandatory to draw a definite conclusion on the nature of 



the transition at physical quark masses Ill4l3l]. 



Among several methods to study the nature of phase transitions, a probability distribution 
of physical observables provides us with the most intuitive way to determine the order of the 
phase transition. Adopting an appropriate physical quantity such as a quark number, chiral order 
parameter, gauge action etc., the corresponding probability distribution function is constructed by 
measuring a generation rate of configurations at each fixed value of the physical quantity. Then 
an existence of double or multiple peaks in the probability distribution function give a signal of a 
first order phase transition, since two phases coexist at a first order phase transition point. 

In this study, we investigate the quark mass dependence of the order of the QCD phase transition 
in the large quark mass region using the hopping parameter expansion. We expect the first order 
phase transition in the heavy quark limit becomes a crossover as we decrease the quark masses 
from infinity. We take the plaquette, i.e. lxl Wilson loop, as the quantity to study the order of 
the transition. The reweighting technique bj, [15|] is employed to vary quark masses in the lowest 
order of the hopping parameter expansion, and we determine an end point where the first order 
phase transition turns into a crossover. 

This paper is organized as follows. Basic properties of the plaquette distribution function 
are discussed in Sec. HU A method to calculate the plaquette distribution function by the hopping 
parameter expansion is introduced in Sec. lIIIl Details of our simulations and results for the effective 
potential defined from the distribution function are presented in Sec. IIVI Results show that the 
first order phase transition becomes a crossover as the quark mass decreases from infinity. We then 
determine the location of the end point for the cases of Nf = 1, 2, 3 and 2 + 1 with the Wilson 
quark action. Our conclusions are given in Sec. IVII 



II. PROBABILITY DISTRIBUTION FUNCTION 

A probability distribution function provides us with one of the most powerful methods to 
determine the order of phase transitions. Since there exist two phases simultaneously at a first 
order phase transition point, we expect that the probability distribution function has multiple 
peaks near the transition point. In this study, we consider the distribution function of the average 
plaquette P, 

(1) 
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where U n „ is the link variable and iV s it e = N$ x Nt is the lattice volume. We adopt the standard 
one-plaquette action (S g ) for glues and the standard Wilson fermion action (S q ) for quarks: 

S g = -6iV site /3P, (2) 
S « = E f " "/ E $° [C 1 " >)^,^S A + (1 + 7m)c1-Am^a1 1 ^ 




(4) 

where iVf is the number of flavors, is the hopping parameter for the /-th flavor, and (3 = 6/g 2 
is the (inverse) gauge coupling. The hopping parameter Kf controls the quark mass, which is 
proportional to for small Kf, while the lattice spacing is mainly controlled by f3. For the case 
of degenerate quark masses, i.e. Kf = K for / = 1, . . . ,Nf, the plaquette distribution function is 
defined as 

w(P',(3,k) = J VUV^jV^ 5(P' - P) e -^-s 9 = I T>U 5{P' - P) (det M{k)) Ni e Gm ^ p . (5) 

The partition function is given by Z(k,(3) = J w{P' , /?, K)dP' . In the followings, we denote the 
argument P' in uu simply as P. 

The plaquette distribution function has the following useful property [151 ]. Under the change 
from /3q to f3, w(P, f3, k) transforms as 

w(P, 0, k) = e^- MN ^ p w{P, /3 , k). (6) 

Therefore, the effective potential defined by 

V eS (P,P,K) = -\nw(P,P,K) (7) 

transforms as 

V cS (P, 0, k) = V eS (P, /3 ,k) - 6(/3 - l3 )N site P. (8) 
From this property, we find that 

^(P,P,k) = ^f(P,f3 ,K) -6(/3-/3 )iV site , (9) 

and d 2 V e fi/dP 2 is independent of /3. 

Since the distribution function is doubly peaked at and around a first order transition point, the 
corresponding effective potential has a double- well structure and the derivative dV e s/dP becomes 
an S-shaped function. Therefore, close to the transition point, dV e s/dP must vanish at three 
points. In general, to find the first order phase transition by observing these properties, a fine 
tuning of /3 is required. For the case of the plaquette effective potential, however, we can detect 
the existence of the first order transition through the appearance of the S-shape in dV e R jdP without 
fine tunings, since the P-dependence of dV e ff/dP itself is independent of /3 as shown in Eq. ([9]). 
After the detection, the first order transition region in /3, in which dV e s/dP vanishes at three values 
of P, can be identified using Eq. (|9|). Therefore, dV e s/dP is useful to determine a region of the 
first order phase transition. 
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We next discuss the K-dependence of V e s considering the ratio of the distribution functions at 
n and kq: 

R(P,K,K ) = r (10) 

fVU5(P' - P)(detM{ K )) N ie^ N ^ p jVU5(P' - P)(det M(k)) n ^ 



JVU5(P' - P)(detM(K )) N !e^ N ^ p JVU5(P' - P)(detM(«o)) 

h(pl _ p , (detM( K )) N ! \ 
°V ^) (det M( K0 ))"t / {/3}Ko) 



(5(P> - P)) 



(11) 



(/3,«o) 



Note that R(P,k,kq) is independent of /3, and thus R(P,k,kq) can be evaluated at any (3. Using 
R(P, k, Kq), the K-dependence of the effective potential is given by 

V cS (P,/3,k) = -lnR(P,K,K ) + V eS (P,f3,Ko). (12) 

This argument can be easily generalized to improved gauge actions including larger Wilson 
loops, for which we redefine the average plaquette as P = —S g /(6N s n e /3). On the other hand, 
/3-dependent improved quark actions make the analysis more complicated. 

We note that the identification of the order of the phase transition by the V e g is equivalent to that 
by the fourth order Binder cumulant, = ((X — (X)) 4 ) / ((X — (X)) 2 } with X an appropriate 
operator signaling the phase transition, since both V e s and B4 detect a change of the distribution 
function. On the other hand, to correctly evaluate B4 at a first order transition point, we need 



to know precisely the probability distribution in a wide range of X covering both peaks 15|. A 
fine-tuning to the transition point and a high statistics with sufficiently many flip-flops are required 
for a reliable estimate of B4 at the transition point. Furthermore, to achieve well separated peaks, 
a large system volume is required. This makes the whole study quite demanding, in particular for 
weak first order transitions such as the case of the heavy quark limit of QCD. On the other hand, 
a reliable V e s in a wide range of P can be easily obtained by combining data at different (3 points 
thanks to the relation ([8]). 



III. QUARK DETERMINANT IN THE HEAVY QUARK REGION 



To investigate the quark mass dependence of the plaquette effective potential, we evaluate the 
quark determinant by a Taylor expansion with respect to the hopping parameter k in the vicinity 
of the simulation point kq: 



In 



det M(k) 
detM(Ko) 



00 1 

T- 

71=1 



9™ (In det M) 



no 



where 



V r , 



d n In det M 

9k" 



(-l) n+1 (n- 1)! tr 



M 



n=l 



(13) 



(14) 



Calculating the derivative of the quark determinant P n , the K-dependence of the effective potential 
can be estimated. 

Adopting kq = 0, the equation (fl4l) reads 



V r . 



l) n+1 (n- 1)! tr 



dM 
8k 



(15) 



K=0 
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/3 configurations bin size 



5.6800 


100,000 


100 


5.6850 


430,000 


2,150 


5.6900 


500,000 


2,000 


5.6925 


670,000 


3,350 


5.7000 


100,000 


500 



TABLE I: The number of configurations and the bin size for jackknife error analyses. 



where M X y = 5 X)V at k = 0. and (dM/dn) Xj y is the gauge connection between x and y. Therefore, 
the nonvanishing contributions to T> n are given by Wilson loops or Polyakov loops. Considering 
the anti-periodic boundary condition and gamma matrices in the hopping terms, the leading order 
contributions to the Taylor expansion are given by 



In 



det M(k) 



detM(0) 

where f2 is the Polyakov loop defined by 



2887V site K 4 P + 12 x 2 Nt N^K Nt ReQ + 



n 



i 



v ^ 1 r 

/ , q tr "n,4^ n+ 4 ) 4^' I1+ 24 ) 4 ' ' - ^n+(JV t -1)4,4 



(16) 



(17) 



The ratio R(P, k, 0) = w(P, (3, k)/w(P, /3, 0) is then calculated at arbitrary values of j3 but small k 
as follows. 



R(P',k,0) 



(6(P' - P) exp[iY f (288A^ site K 4 P + 12 x 2 Nt N^n Nt Ren + •••)]> 



(/9,«o=0) 



(5(P> - P)) 



(/3,ko=0) 



= 288Ar f Af sito K 4 P' 



, (5{P' - P) exp[JV f (12 x 2 Nt N^n Nt ReQ + •••)]> 



03ao=o) 



(5(P> - P)) 



(/3,« =0) 



(18) 



For the case of Nt = 4 we study, the truncation error is (9(k 6 ). The contribution from the plaquette 
can be absorbed by a shift of /3. 



IV. NUMERICAL SIMULATIONS AND THE RESULTS 
A. Simulation parameters and plaquette effective potential 

In the heavy quark limit we perform simulations of SU(3) pure gauge theory on a 24 3 x 4 lattice. 
We generate the configurations by a pseudo heat bath algorithm followed by four over-relaxation 
sweeps. We perform simulations at five (3 values in the range 5.68 - 5.70. Simulation points and 
statistics are summarized in Table HI 

The results for the histogram of P, i.e. the plaquette distribution function Eq. ([5]) normalized 
by the partition function Z, are plotted in Fig. [TJ In this calculation, we approximate the delta 
function by a Gaussian function: 5(x) « l/(A- v /7r) exp [— (x/A) 2 ], where A corresponds to the 
width of the Gaussian function. As A decreases, the resolution of the distribution function becomes 
better, while the statistical error increases. Examining the resolution and the statistical error, we 
adopt A = 0.000283. This figure shows that the value of P generated in a simulation at single f3 
distributes in a narrow range. 
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FIG. 1: Plaquette histogram w/Z obtained by 
SU(3) pure gauge simulations at f) = 5.68-5.70 
on the 24 3 x 4 lattice. 
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FIG. 2: Derivative of the effective potential in 
SU(3) pure gauge theory at /3 — 5.69 converting- 
data obtained at j3 = 5.68-5.70 by Eq. ©. The 
bold black curve is an average over the results 
obtained at different j3. 
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FIG. 3: Derivative of the effective potential at 
nonzero k in two-flavor QCD. 
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FIG. 4: Effective potential V c g at /?trans for each k. 



We then calculate the derivative dV e g/dP by the difference between the potentials at P — ep/2 
and P + ep/2. The value of ep is set to be 0.0001, considering the resolution in P. The ep- 
dependence is much smaller than the statistical error around this value of ep. Results of dV c f[/dP 
at k = are shown in Fig. [21 where data points with low statistics are removed for clarity. In 
this figure, we adjust results at different /3's to (3 = 5.69 by using Eq. ([9]). Results of dV e ^/dP 
obtained in simulations at different f3 are consistent with each other within statistical errors, though 
the ranges of P in which V e s is reliably obtained are different. The jackknife method is used to 
estimate the statistical error of the effective potential and its derivatives. The bin size of the 
jackknife analysis is listed in Table HI We then combine these data by taking an average weighted 
with the inverse-square of errors in the overlapping regions. We here exclude data at P far away 
from the peak of the distribution at each f3 which thus has poor statistics. The final result for 
dV c g I dP at k = is shown by the black line without symbols in Fig. [2l We find that dV e s / dP 
is not a monotonically increasing function at n = 0, indicating that the effective potential has a 
double- well shape. 

The quark mass dependence of the effective potential is investigated by calculating R(P, k, 0) 
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FIG. 5: Curvature of the effective potential at k = 0.058 (left), 0.062 (middle), 0.066 (right). Black symbols 
are results at K = 0. Bold horizontal lines on the line of d 2 V e g jdP 2 = represent errors of the plaquette 
values where d 2 V c e/dP 2 vanish. 



up to the order k 4 in Eq. (|18p . Using the data of w(P,(3,k = 0) and R(P,k,0), we evaluate the 
dV c s/dP at nonzero k. Results for Nf = 2 are plotted in Fig. [3j The S-shape structure becomes 
weaker as k increases and turns into a monotonically increasing function around k = 0.066. This 
behavior suggests that the first order phase transition at k = becomes weaker as k increases 
and the transition becomes a crossover at k ~ 0.066. In Fig. [J] we show the k dependence of V c r 
obtained by a numerical integration of dV c s/dP. In this figure, (3 is adjusted such that the two 
minima have the same height (see Sec. II V B 21 for details), and the integration is started from the 
peak point P pea k of V e s between the two minima. We find that the double-well becomes shallower 
as k increases and become almost flat around n = 0.066. 



B. Critical point in the heavy quark region for Nf = 2 

In this subsection, we evaluate the value of n at which the first order phase transition terminates. 
We denote the corresponding critical point as n cp . 



1. Critical point from d 2 V c s / dP 2 

We first calculate the second derivative of V c g(P, /?, n) by a numerical differentiation of dV c g/dP. 
When the transition is of first order, there is a region in P where d 2 V e g / dP 2 is negative between 
the two bottoms of V e g (P, {3, k). The first order transition region is thus identified by calculating 
the sign of the second derivative of V^ff- As discussed in Sec. HH the second derivative of V e s(P, (3, k) 
is independent of (3. Therefore, the identification can be performed at any f3. 

In Fig. we plot the results of d 2 V eS /dP 2 at k = 0.058 (left), 0.062 (middle) and 0.066 (right) 
for Nf = 2, together with the results at k = which are shown as black symbols. We find that 
the region where d 2 V e g/dP 2 < becomes narrower as n increases. Bold horizontal lines at zero 
show the range of P where the curvature vanishes within statistical errors. Hereafter we denote the 
position at which d 2 V e g/dP 2 = as Pjy" =0 )(k) or K(y» =0 )(P). Results of -P(y'^=o) are plotted in 
Fig. [6] as a function of k. The region of the negative curvature seems to vanish around k = 0.066. 
Since the curvature is always positive beyond this k, the effective potential is no longer a double- well 
type, i.e. the transition becomes a crossover at k ^ 0.066. 

To evaluate the critical point n cp we study the K-dependence of d 2 V e g/dP 2 at fixed P, as shown 
in Fig. [7] for typical points. Data with large statistical errors are removed from this figure. Fitting 
data by a linear function in k, we obtain K^y^ =0 ^p-j, which is shown by a square symbol in the 
same figure. The horizontal bar represents the statistical error. Results of K(y^ =0 ) are plotted in 
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FIG. 6: Value of P where d 2 V c g/dP 2 vanishes at each k. 
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FIG. 7: K-dependence of d 2 V e g/dP 2 and the results 
of K ( y c » =0 ) at P = 0.5479, 0.5483 and 0.5486. 
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FIG. 8: Value of k where d 2 V c s/dP 2 vanishes at 
each P. 



Fig. [8]for 0.5474 < P < 0.5487. When k is larger than the maximum value of ^(v" a =o)i d 2 V e s/dP 2 
is non- negative for all values of P in this region. Therefore, the maximum value of H(y" s =o) i s 
nothing but the critical point K cp . We obtain K cp = 0.0685(72) for TVf = 2. 

2. Critical point from the double-well structure of V c g 

Alternatively, we may determine K cp as the point where the two minima of V c s(P) merge and 

the barrier between the two minima vanishes. The double- well structure of V e s is most clearly 

seen at (3 where the two minima of V c s has the same height (see Fig. d]). Such f3, say /3trans> can 

be determined by a Maxwell's construction for dV e s/dP: Let us denote the values of P at the two 

l ' B dV cS 

-{P,p trans ,K) dP 

'Pa 

or, equivalently, 



f 1 

minima as Pa and Pb- The condition V c s(Pa) = V c s(Pb) implies / 

Jp, 



dP 



0, 



.PA 



P° ak dVeS 

dP 



(P,P trans ,K)dP 



Pb 



peak 



dVoS 
dP 



(P,^ ns ,K)dP 



(19) 



where -P pea k is the peak position of V e g between Pa and Pb at which dV e g/dP vanishes. Changing 
(3 by Eq. ([9]), we search /3 trans where Eq. ([19]) is satisfied. The results for /3 trans are plotted in the left 
panel of Fig. [9] as a function of k. At n = 0, we obtain /3trans = 5.69138(3). At ft = 0.066, because 
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FIG. 10: Potential barrier height Vp ea k (Left) and the gap AP (Right) between the two minima at ptrans 
for Nf = 2. Critical point n cp is estimated by linear extrapolations shown in the figures. 



the minima of V e g are too shallow to apply the Maxwell construction, we instead determine /3t ra ns 
as the point where the bottom region of V e s becomes flat, as shown in Fig. 01 In the right pane of 
Fig. [9l we replot Ptrans as a function of k 4 , as motivated by the hopping parameter expansion. We 
note that the data is well fitted by a linear function of k . 

In an approximation to disregard fluctuations around the plaquette expectation value P, we 
can identify V e s(P) with the free energy of the system. Then the condition V e s(P\) = V c s(Pb) 
means that the pressure p is balanced between the phases A and B. Therefore, we may view /3t ra ns 
as an estimate of the first order transition point. Actually, our result /3t ra ns = 5.69138(3) for k = 
is quite close to the first order transition point /?trans = 5.69153(86)-5.69211(35) defined by the 
peak position of the plaquette susceptibility in SU(3) pure gauge theory with the same lattice size 
24 3 x 4 0. 

In Fig. [101 results of the potential barrier height Vpeak = Kff(-Ppeak) — ^eff(-f*A) and the gap 
AP = Pb — Pa at f3 trans are shown as functions of k. We find that both V pea k and AP decrease 
drastically above k = 0.05 and vanish around k ~ 0.066. Performing a linear extrapolation using 
three nearest points of Vp ea k [AP] as shown in Fig. [10J we obtain K C p — 0.0647(6) [0.0662(4)]. A 
corresponding value of /3 cp at the critical point can be estimated by extrapolating /3t rans to fc C p, 
assuming a linear function in k or k 4 as plotted in Fig. [9j Results for K cp and /3 cp are summarized 
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in Table ITTl In the table, we also present /3 cp for K cp determined by d 2 V e g/dP 2 . From these results, 
we obtain 

K cp = 0.0658(3) (j^), p cp = 5.6819(1) (5) (20) 

where the central values and their statistical errors given in the first brackets are determined by a 
weighted average of the results, and systematic errors given in the second brackets are determined 
from the scattering of the results due to the method and extrapolation function, neglecting the 
data from the d 2 V c g/dP 2 method which are consistent with other results within large errors. 

To get a rough idea about the value of the pseudoscalar meson mass (mps) corresponding to the 
critical point, we estimate the ratio T/mps = l/(N t mpsa) at the critical point. Using a data of 
zero-temperature pseudoscalar meson mass in Nr = 2 QCD along the finite temperature crossover 
curve for Nt = 4 in the range k = 0.16-0.19 [161], we perform a fit (mpso) -1 = f\K 2 + /2^ 4 with 
fitting parameters f\ and /2, where the constraint (mpso)" 1 = at k = is taken into account. 
We find T/mps ~ 0.023 for Nf = 2 at the critical point n cp « 0.066. Because T trans /mps = 0(1) 
around the physical point, the small value 0.023 means that the critical pseudoscalar meson mass 
is much larger than the physical pion mass. 







/3 CP 


method 




k fit k 4 fit 


^pcak 


0.0647(06) 


5.6824(02) 5.6823(03) 


AP 


0.0662(04) 


5.6818(01) 5.6814(02) 


d 2 V eS /dP 2 


0.0685(72) 


5.6808(30) 5.6798(50) 


total 


0.0658(03) (tu) 


5.6819(1)(5) 



TABLE II: Critical point K cp and /3 cp defined by Vp Ca k, AP and d 2 V e g/dP 2 



C. Latent heat 

The gap of the internal energy density, Ae, at a first order transition point is called the latent 
heat. Since p is continuous there, Ae can be evaluated as the gap of the trace anomaly: 

e — 3p 1 dlnZ rm 

[T = 0] 



T 4 VT 3 da 



K 



a^-6(P) + iV f o^ (1152k 3 (P) + AT 1 12 x 2 Nt z^ 1 (ReO)) + 0( 
da da ' ' 

- 0] (21) 



where [T = 0] is the zero temperature contribution to be subtracted for renormalization. At k = 0, 
Qj6)a(df3/da) « -(0.064-0.078) at j3 = (3 c (N t = 4), depending on the method to define the scale 



17h2Q|]. At k = 0, (Refi) = due to the center Z{3) symmetry, and a(dn/da) = because k = 
is a line of constant physics. Neglecting terms proportional to K 3 a(dn/da) at small k, the latent 
heat is roughly evaluated as 

^ = A(£ - 3P) ^-6A^(AP) (22) 
T 4 T 4 1 da x ' v ; 

where AP is shown in Fig. [10] for Nf = 2. It is found from the behavior of AP that the latent heat 
decreases as k increases. 
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0.05 0.1 



Ni 




1 

2 
3 


0.0783(4) (t$ 3 ) 
0.0658(3) (±f x ) 
0.0595(3) (±* ) 



TABLE III: K cp for N t = 1, 2 and 3 deter- 
mined by the reweighting method with the 
lowest order hopping parameter expansion at 
At = 4. 



FIG. 11: The phase boundary separating the 
first order transition region and crossover re- 
gion in the (re u dj K s) plane. 



V. THE CASES OF 2 + 1-FLAVOR AND DEGENERATE A f -FLAVOR QCD 

The analysis can be easily generalized to an arbitrary value of Nf and also to the case of 
Nf = 2 + 1 QCD. At the leading order of the hopping parameter expansion, the contribution of 
quark determinants in the partition function is given by 



In 



(detM(K ud )) 2 detM(/« s ) 



(detM(O)) 3 



288i\U,(2«4j + 4)P + 12 x 2 Ar 'iV s 3 (2^ + /e^)Refi + • • • (23) 



for the case of iVf = 2 + 1, where K ud and k s are hopping parameters for light and strange quarks. 
The modification factor for the reweighting in k thus becomes 



R(P',K nd ,K s ,0) = e 288iV 8i te(2< d +<)P' 

5{P' - P) exp[12 x 2^JV a 3 (2«^ + ^)ReQ + ■■■}) 



(24) 



(/3,«=0) 



Since the contribution from the plaquette in this equation does not affect the second derivative of 
V e fi, the difference from the case of Nf = 2 is just the replacement of 2n Nt by 2k^ + k^*. Thus, 
the line which separates the first order phase transition and the crossover is given by 

2(« lld ) J * + = 2«'= 2 )^ (25) 

where N t = 4 and k^= 2 = 0.0658(3) (1^) as determined in Sec. HVBl We draw the line in Fig. [Til 
in which the colored region corresponds to the first order transition. The case of degenerate N{- 
flavor QCD can be similarly investigated. Results for K cp are summarized in Table ITTTI for Nf = 1, 2 



and 3. Our K cp for Nf = 1 is consistent with the result obtained by an effective model in Ref. 2ll j. 
In the present approximation of the lowest order hopping parameter expansion at Nt = 4, /3 cp is 
independent of Nf. 
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VI. CONCLUSION 



We have investigated the order of the deconfinement phase transition in the heavy quark region 
of QCD by calculating an effective potential V e s(P) defined as the logarithms of the probability 
distribution function for the average plaquette P. To study the fate of the first order deconfinement 
transition, V^ff has to be reliably evaluated in a wide range of P covering the both phases. Applying 
a reweighting method, we combine the results at five different /3 points to calculate V e s in a wide 
range of P. We evaluate V e g at k = by simulations in SU(3) pure gauge theory, and reweight 
it to nonzero but small values of k using the leading order hopping parameter expansion on an 
N t = 4 lattice. 

At k = 0, V^ff(P) show clear double- well structure in accordance with the first order deconfine- 
ment transition of the SU(3) pure gauge theory. When we increase k from zero, the double-well 
shape becomes weaker, and eventually disappear at finite k, indicating that the first order transi- 
tion turns into a crossover at that point as suggested from the effective Z(3) model with an external 
magnetic field. We estimated the critical point K cp by examining various properties of V e f[. The 
results for the critical point in degenerate Nf -flavor QCD as well as that in the Nf = 2 + f QCD 
are summarized in Table IIIII and Fig. [IT] 

The calculations are done in the leading order approximation of the hopping parameter expan- 
sion on an Nt = 4 lattice. Although the values of K cp we obtained are quite small, it is important 
to confirm the reliability of the leading order approximation quantitatively. In the next leading 
order, we have 0(k 6 ) contributions from Wilson loops with length 6 and 0(n Nt+2 ) contributions 
from Polyakov loops which have N + 2 link variables including two spatial links. To estimate 
the effects of them, measurement of the probability distribution function including expectation 
values for these operators is required. Another important point to be studied is the continuum 
extrapolation by increasing Nt- We leave these studies as future works. 
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